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Abstract. Software implementations of controllers for physical systems are at the core of many embedded 
systems. The design of controllers uses the theory of dynamical systems to construct a mathematical control 
law that ensures that the controlled system has certain properties, such as asymptotic convergence to an 
equilibrium point, while optimizing some performance criteria. However, owing to quantization errors arising 
from the use of fixed-point arithmetic, the implementation of this control law can only guarantee practical 
stability: under the actions of the implementation, the trajectories of the controlled system converge to a 
bounded set around the equilibrium point, and the size of the bounded set is proportional to the error in the 
implementation. The problem of verifying whether a controller implementation achieves practical stability 
for a given bounded set has been studied before. In this paper, we change the emphasis from verification to 
automatic synthesis. Using synthesis, the need for formal verification can be considerably reduced thereby 
reducing the design time as well as design cost of embedded control software. 

We give a methodology and a tool to synthesize embedded control software that is Pareto optimal w.r.t. 
both performance criteria and practical stability regions. Our technique is a combination of static analysis to 
estimate quantization errors for specific controller implementations and stochastic local search over the space 
of possible controllers using particle swarm optimization. The effectiveness of our technique is illustrated using 
examples of various standard control systems: in most examples, we achieve controllers with close LQR-LQG 
performance but with implementation errors, hence regions of practical stability, several times as small. 



1. Introduction 

Software implementations of controllers for physical systems are the core of many critical cyber-physical 
systems. The design of these systems usually proceeds in two steps. First, starting with a mathematical 
model of the system, one designs a mathematical control law that ensures that the physical system, equipped 
with this control law, has certain desirable properties such as asymptotic stability (convergence to an ideal 
behavior) and performance. Second, the control law is implemented as a software task on a specific hardware 
architecture. Since the implementation has quantization errors due to the use of fixed-precision representation 
of real numbers, the quantization of a stabilizing controller may lead to limit cycles and chaotic behavior 
[Kal56| . Hence, the implemented system usually guarantees the weaker property of practical stability, where 
the system is guaranteed to converge to a bounded set around the ideal behavior and the size of the bounded 
set is proportional to the quantization error. 

Much recent research has focused on verifying that a given implementation of a control law guarantees that 
the practical stabihty region hes within a given set |PW06[ iFWFZl IFerlOi lAMSTlOi fDMllI . In this paper, we 
change the emphasis from verification to synthesis. We provide a design methodology to synthesize a control 
implementation for which the effect of implementation errors on system performance is minimized. Hence, the 
need for verification can be substantially reduced. 

We focus on linear systems in this paper. For linear systems, a standard optimal control design approach uses 
the linear quadratic regulator (LQR) and linear quadratic Gaussian (LQG) algorithms [Hcs09 , which finds a 
feedback controller stabilizing the plant while minimizing quadratic cost functions. The LQR cost function 
takes into account the deviations of the state and control inputs from ideal values and the LQG cost function 
takes into account the deviation of the state from its estimation. However, in general, they do not take 
implementation errors arising from fixed-precision arithmetic into account. Thus, a controller optimizing only 
the LQR-LQG cost may have a large implementation error because its implementation on a fixed-precision 
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platform has large numerical errors, but a controller "close" to the optimal performance may have much lower 
numerical errors when implemented on the same platform. 

In our methodology, we modify the LQR-LQG performance criterion to additionally minimize the error due to 
quantization in the implementation. Technically, our methodology has two parts. First, how can we estimate 
the quantization error of a given implementation? Second, how can we find Pareto-optimal points for the two 
objectives given by the LQR-LQG and quantization error cost functions? We proceed as follows. 

For the first step, for a given linear feedback controller and the operating intervals of the states of the plant 
and the controller, we first perform a precise range analysis of the controller variables, and use the computed 
ranges to allocate bitwidths to each controller variable. We implement our range analysis based on linear 
programming. Using the allocated bitwidths, we generate code for a fixed-precision program implementing 
the control law. Finally, we use an algorithm based on mixed-integer linear programming to find a bound on 
the maximum difference between the ideal control law and the output of the fixed-precision program. 

For the second step, we optimize a weighted linear combination of the two cost functions using a stochastic 
local search technique. LQR-LQG is attractive because it gives rise to a convex optimization problem, for 
which efficient solutions are known. Unfortunately, additionally tracking the quantization error results in a 
non-convex optimization problem. We solve the non-convex optimization problem using particle swarm opti- 
mization (PSO), a population-based stochastic optimization approach |KE95[ |LAS09[ ITLY07| . PSO iteratively 
solves an optimization problem by maintaining a population (or swarm) of candidate controllers, called par- 
ticles, and moving them around in the search-space of possible controllers, trying to minimize the objective 
function. 

In more detail, our algorithm proceeds as follows. Given a linear control design problem, we set up a non- 
convex optimization problem to minimize a weighted combination of the LQR-LQG cost function and the 
implementation error. We minimize this cost function using PSO. In each step of PSO, given a new position 
of a particle, we check if the position represents a stabilizing controller (by examining the eigenvalues of the 
controlled system). If not, we assign the position an infinite cost. If the position represents a stabilizing 
controller, we generate the best possible fixed-point code for this controller under a hardware budget and 
perform static analysis to estimate a bound on the implementation error. We compute the value of the 
objective function by taking the weighted sum of the LQR-LQG cost and this bound. We continue PSO until 
convergence or until some iteration bound is met. At this point, we output the controller that minimized the 
objective function. 

We have implemented this methodology on top of Matlab's Control Theory Toolbox, using an implementa- 
tion of PSO proposed in [EKGJL2 , and a custom static analysis using the lp_solve linear programming tool. 
In our experiments, we compare the LQR-LQG cost and implementation errors of controllers generated by 
conventional LQR-LQG optimization (implemented in Matlab) with controllers generated by PSO using our 
methodology. In most cases, our controllers have LQR-LQG costs close to the optimal LQR-LQG controllers, 
but have implementation errors that are reduced by a factor of 4 or more. Thus, we generate controllers with 
guaranteed bounds on practical stability regions that are 4 times or more smaller than the pure LQR-LQG 
controllers. Our work provides, for the first time, an integrated analysis and tool to take quantization errors 
into account in model-based design and implementation of controllers. 



Other Related Work Besides the related work mentioned above, we mention the results in |Wil85[ IWil891 
ILSG 92J which provide controller synthesis approaches minimizing some performance criteria while controllers 
are implemented using fixed-point arithmetic. The results in |Wil85| IWil89| ILSG92| assume some excitation 
conditions under which the quantization error can be modeled as a zero mean uniform white noise. Further- 
more, they do not provide any bounds on regions of practical stability. However, the result in this paper 
does not have any assumption on the quantization error and it provides an explicit bound on the regions of 
practical stability. 
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The range analysis problem has been studied extensively in the context of optimum bitwidth allocation to 
intermediate variables in a fixed-point program, mostly in the DSP domain. Both static |LGC+06[ ILCNT07I 
IOCC+OT] and simulation-based }BR05[ IMSBZ07] approaches have been used. Static approaches usually em- 
ploy abstractions based on interval arithmetic |Moo66| or affine arithmetic |SF97] . Simulation-based methods, 
especially those performing constrained-random simulations, suffer from the lack of completeness. This causes 
the resultant system to be non-robust, incomplete simulation can lead to overflow conditions resulting in in- 
correct behavior. We found our mixed-integer linear programming approach to be both precise and reasonably 
fast for our application. 



2. Preliminaries 



2.1. Controllers and Observers. We use symbols Nq, K, and Rq for the set of nonnegative integers, real 
and nonnegative real numbers. For a vector x £ ]R", we denote by Xi the i-th element of x, and by the 
Euclidean norm of x. Recall that = + + • • • + a;^ . The symbols /„ and Onxm denote the identity 
and zero matrices in M"^" and M"^™, respectively. A continuous function 7 : M.q — > M.q , is said to belong 
to class /C if it is strictly increasing and 7(0) = 0; 7 is said to belong to class /Coo if 7 G /C and j{r) — > cxd as 
r — 00. A continuous function /3 : Rq x R'^ R^ is said to belong to class JCC if, for each fixed s, the map 
/3(r, s) belongs to class /Coo with respect to r and, for each fixed nonzero r, the map /3(r, s) is decreasing with 
respect to s and /3(r, s) — > as s — > 00. 

In this paper, we focus on linear control systems given by the differential equation: 

^ = A^ + Bv + Buj, 



(2.1) 



where, for any t G K, ^(t) G R", v{t) G E™, u{t) G K«, r}{t) G R^, and A, B, B, and C are matrices of 
appropriate dimensions. The curve ^ : E — ?> E" is a trajectory of (2.1) if there exist curves u : E — ^ E™ and 
w : E — >■ E'^ such that the time derivative of ^ satisfies (2.1). In the rest of the paper, we assume that all 
curves v and uj have some regularity assumptions, guaranteeing existence and uniqueness of the solutions of 
(2.1). Note that v, w, 77, and v denote control input, disturbance, output of the system and measurement 
noise, respectively. We assume that uj{t) and v{t), for any i G E, are zero-mean Gaussian noise processes 
(uncorrelated from each other). For all curves w, we also write S,xv{t) to denote the points reached at time t 
under the input v from initial condition x = ^xv{^)- 

To describe the mismatch between the controller specifications and its software implementations such as digital 
sampling and finite precision arithmetic, which is the focus of this paper, we consider the discrete-time version 
of (2.1 1, as follows: 



(2.2) 



x[r + 1] = Arx[r] -I- Bru[r] + Brd[r] 
y[r] = Cx[r] -\- v[r], 



where the matrices A^., B^., and B-,. are given by: 

r.(r-hl)r /-(r+ljT 

Aq- ^ e , -S-r ~— 



/ e^^^-'^Bdt, e^^^-'^Bdt, 

J rr J rr 



and T is the sampling time. The function e"**, for any t G E, denotes the matrix function defined by the 
convergent series: 

e^*=/„+At+iA2i2 + i.A3t3 + ... , 

where e is Euler's constant. The signals x, u, d, y, and v describe the exact value of the signals ^, v, lo, ry, and 
V, respectively, at the sampling instants 0, t, 2r, 3t, .... Mathematically, we have: 

x[r\ — £,{rT), u[r] = v{rT), d[r] ~ uj^rr), y[r] — r]{rT), v[r] — ly^rr), Vr G Nq. 



The term Cs in (2.2) is the sampling error. It can be shown that by sampling sufficiently fast, the error 
can be made arbitrarily small ICF95j . Since typical embedded controller implementations use sampling time 
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in the range of milliseconds to microseconds, we will make the assumption that quantization errors dominate 
the sampling errors, and assume that — 0. 

We assume that only output of the system y is measurable and not the full state x. Hence, a (proportional) 
feedback K : M" — >■ K™ defines the input u[r] = —Kx[r] based on an estimation x of the state x. As explained 
in |Hes09| . the estimation x can be constructed using the observer dynamic: 



(2.3) 



x[r + 1] = ^ra;[r] + Bru[r] + L {y[r] — y[r]) , 
y[r] = Cx[r], 



where y should be viewed as an estimate of y and the linear map L : — ?> is called an observer gain. By 
applying the feedback u[r] = —Kx[r] and combining the dynamics of control system in (2.2) and observer in 



(2.3), one obtains: 



(2.4) 



x[r + 1] = Arx[r] - BrKx[r] + Brd[r], 

x[r + 1] = {Ar - BrK - LC)x[r] + LCx[r] + Lv[r]. 



As shown in |AMST10] . using a fixed-point implementation of the feedback gain as well as the observer 
dynamic, one gets the following overall dynamic: 

x[r + 1] = ^ra;[r] - BrKx[r\ + Brd[r] + Breq2, 
x[r + 1] = {Ar - BrK ~ LC)x\r\ + LCx\t\ + Lv\r\ + e,i. 



(2.5) 



where e^i and 6^2 are quantization errors in observer dynamic and feedback gain codes, respectively. Now, 
one can rewrite the control system in (2.5) as follows: 



(2.6) 
with: 

and: 



w{r + 1] = Gw{r\ + i?iei[r] + i72e2H, 





X 




■ d ' 




Gql 


w = 




, ei = 




, 62 = 






X 




V 


6,2 



G = 



Ar 



BrK 



LC Ar - BrK - LC 



Hi = 



Br 

Onx g 



L 



Ho = 



-'nx n 



Br 



Since states of the control system (2.1 ) are bounded physical quantities, such as temperature, pressure, position, 
velocity, acceleration and so on, their estimations and the output of the control system are bounded quantities 
as well. Hence, in the rest of the paper and without loss of generality, we assume that y e y, and x & X, 
where Y C M^, and X C M" are compact. 



2.2. Stability of perturbed systems. Here, we recall the notion of uniform global asymptotic stability with 
respect to a set, presented in |LSW96| . 



Definition 2.1 f |LSW96j ). A control system of the form (2.1) is uniformly globally asymptotically stable 



(UGAS) with respect to a set A if there exists a ICC function f3 such that for any t S Mq , any x € 



any 

D2 C K'J, where Di, and D2 are 



control input v : Mq — > Di C R™, and for all possible disturbances uj : 
compact sets, the following condition is satisfied: 

(2.7) u..mA<P{M\A,t), 

where the point-to-set distance ||a;|l_4 is defined by \\x\\j( = infj,g_4 ||x — y\\. 

When the set ^ is a singleton {xq}, we speak of an asymptotically stable equilibrium point xq rather than a 
UGAS set. The notion of UGAS for discrete-time control systems is obtained from Definition [O] by replacing 
t e Mfj" with r € Nq. 



We recall the following result describing how stability properties are affected by additive disturbances. 
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Proposition 2.2 ( [AMSTIO] ). Consider the discrete-time linear system: 

x[r + 1] — Ax[r], 

and assume that the origin is an asymptotically stable equilibrium point. Then, for any signal d : Nq 
satisfying \\d[r]\\ < b{d) for any r G No and some constant b{d) G , the system: 



(2.8) 

is UGAS with respect to the set: 

where 7 is given by: 
(2.9) 



x[r + 1] = Ax[r] + Bd[r], 
^ = {.TGM" I 11x11 < 76(d)}: 



7 = max 
ee[o, 2ir[ 



[e^'l,,-A)-'B 



with i = \/— T. Moreover, the output y — Cx is guaranteed to converge to the set: 

(2.10) Ay^{yeW \\\y\\<lybid)}, 

with: 



"fy = max 
ee[o, 2Tr[ 



C{e''ln-A) 



In control theory, 7j, is known as the £2 gain of the control system in (2.8) with the output y — Cx. The 
following proposition follows from Proposition 2.2 and describes the stability properties of linear control 



systems in (2.6) with respect to disturbance, measurement noise, and implementation errors in the feedback 



gain and observer dynamic. 



Proposition 2.3. Consider the discrete-time linear system in {2.6). Then for any input ei and 62 satisfying 
llei[?']|l < b{ei) and |le2[r]|] < 6(62) for any r G No and some constants 6(ei),fe(e2) G MJ, the system is UGAS 
with respect to the set: 

A = {xe M" I 11x11 < 7i6(ei) + 72^62)} , 

where 71 and 72 are given by: 



7o = max 
ee[o. 2tt[ 



, for 1,2, 



with i = \J —\. Moreover, the output W 9 y = [C Opxn] w is guaranteed to converge to the set: 
(2.11) Ay^{y^W\ \\y\\ < ^^yb{e,) + 72^6(62)} , 

where "/ly and "f2y given by: 



(2.12) 



7,u = max 



[C Opx„] (e^'hu - G) 



, for j = 1,2. 



The error vector ei includes disturbance and measurement noise, depending for example on the environment 
and the quality of the sensors collecting measurements. Hence, the controller designer does not have any 
control on the value of b{ei). However, one can reduce the amount of 71 j, by appropriately choosing gains K 
and L. On the other hand, one can reduce the amount of not only j2y but also 6(62) by appropriately choosing 
gains K and L. We use Proposition |2 .3| in the following way. Given a feedback gain K and an observer gain L, 
we compute £2 gains ^ly and 721, and an upper bound 6(62) on the implementation error 62. Then the output 
of the controlled system (with implementation error) must converge to set Ay in (2.11). We show later that 



appropriate choices of gains K and L can shrink the size of the set Ay and hence, provide a tighter bound on 
the set to which the output of the system converges. 
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2.3. LQR-LQG performance. In addition to asymptotic stability, controller designers also consider the 
performance of the controller, that is, of the controllers ensuring asymptotic stability of the origin, one desires 
the controller that minimizes a given cost function. A common approach for optimal output feedback controller 
are the linear quadratic regulator (LQR) and linear quadratic Gaussian (LQG). The LQR cost function to be 
minimized is given by: 



(2.13) 



Jlqr 



+ CXD 

E 

r=0 



{x[rf Qx[r] + u[rf Ru[r]} 



for some chosen weight matrices Q and R that are positive definite and of appropriate dimensions. 
The LQG cost function to be minimized is given by: 



(2.14) 



Jlqg 



lim E 



where E stands for expected value and e is the estimation error for the control system in (2.4) whose dynamic 
is given by: 

(2.15) e[r + I] ^ x[r + 1] - x[r + 1] ^ {Ar ~ LC)e[r] + Brd[r] - Lv[r]. 

As mentioned before, d and v are assumed to be zero-mean Gaussian noise process (uncorrelated from each 
other) with covariance matrices: 

(2.16) E {d[r]d[rf) = Q, E {v[r]v[rf) = R, Vr e Nq, 
where Q and R are some positive semi-definite matrices of appropriate dimensions. 



A standard control-theoretic construction rewrites the cost function (2.131 as Jlqr = xlO]"^ S{K)x[0], where 
u = —Kx, and S{K) e M"^" is a positive definite matrix that is the unique solution for S to the Lyapunov 
equation: 

(2.17) {Ar - BrKf S {Ar - BrK) -S + Q + K^RK = 0, 

where K is a. controller making Ar — BrK Hurwitzj^ See [Hes09j for detailed information. Additionally, we 
have 



(2.18) 



A„,in(^(K))|la;[0]|p < J, 



LQR 



< A, 



.{S{K))\\x[0]r, 



where Xniin{S{K)) G M+ and Xmax{S{K)) e M+ are minimum and maximum eigenvalues of S{K), respectively. 
Therefore, Jlqr can be minimized for all possible choice of initial conditions by just minimizing the maximum 
eigenvalue of S(K). Note that since 5 is a positive definite matrix, its maximum eigenvalue is equal to its 
induced 2 norirri 11511 . 

Similarly, the cost function (2.14) can be rewritten as Jlqg = ||-P(^)ll7 where P{L) G M"^" is a positive 
definite matrix that is the unique solution for P to the Lyapunov equation: 

(2.19) {Ar - LC) P {Ar - LCf -P + BrQB^ + LRLF = 0, 

where L is an observer gain making Ar — LC Hurwitz. See jHesOQ] for more detailed information. Therefore, 
Jlqg can be minimized by just minimizing ||P(L)||. 



Note that the optimal feedback u — ~Kx minimizing the LQR cost in (2.13) is computed using the determin- 
istic dynamic: 

x[r -|- 1] = A7.a;[r] -f Bru[r]. 



On the other hand, the optimal gain L minimizing the LQG cost in (2.14) is computed using the stochastic 



dynamic in (2.15). Thanks to the separation principle for linear control systems [Hes09| . one concludes that 
the overall closed loop system in (2.4) is UGAS even though the gains K and L are designed separately. 



^ We call the matrix Ar — BrK Hurwitz if its eigenvalues are inside the unit circle, centered at the origin. 
^We recall that induced 2 norm of a matrix A e K"^™ is given by: = ^/Amax (A'^A). 
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Figure 1. Evolution of the output y with initial state (0.2, 0.2)^ for the pair of gains Ki, 
Li and K2^ L2 using 16-bit implementation. 



2.4. The effect of errors. Example We now present a simple motivating example showing how different 
choice of controllers result in different steady state errors due to their fixed-point implementations, yet provid- 
ing approximately the same LQR-LQG performance. Consider a simple physical model of a bicycle, borrowed 
from |AM08j . The dyn amies of the system is given by: 



(2.20) 



6 
6 



"0 f 






■ Ci 


1 








2 " 

ava 






bh bh 




. ^2 . 



v, 



where ^1 is the steering angular velocity, ^2 is the steering angle, 77 is the role angle, v is the torque applied to 
the handle bars, g = 9.8m/s^ is the acceleration due to gravity, h = 1.5m is the height of the center of mass, 
Vq = 2m/ s is the velocity of the bicycle at the rear wheel, a = 0.5m is the distance of the center of mass from 
a vertical line through the contact point of the rear wheel and b = Im is the wheel base. 



The control objective is to design a feedback gain K e 



such that 

the feedback control law u — —Kx, where x— [xi, X2]'^ is the state of the observer in (2.3), makes the 



and an observer gain L G 



closed loop system U GAS. By choosing the matrices Q = I2 and R = 1 inside the LQR cost function and 
Q = 1 and i? = 1 in (2.16), the feedback and observer gains minimizing the LQR and LQG costs are given 
by Ki = [5.1538, 12.9724], and Li = [0.0317, 0.0118]^, respectively. Consider a second pair of feedback 
and observer gains given by K2 = [3.0253, 12.6089] and L2 = [0.0132, 0.1021]^. For the initial condition 
X = (0.2, 0.2)^, the value of the LQR cost function is 264.1908 for feedback gain Ki and 284.1578 for K2. 
Moreover, the value of the LQG cost function is 0.0229 for observer gain Li and 0.0246 for L2- That is, the 
gains K2 and L2 give cost functions about 7% greater than the optimal gains Ki and Li. 

We now show how different choice of feedback and observer gains result in different fixed-point implementation 
errors. For now, let us assume that uj{t) = and i^{t) = 0, for any t S Rj|. In Figure [l] we show the output 
of the closed-loop system starting from the initial condition x = (0.2, 0.2)^, when the feedback gain and 
observer dynamic are implemented using 16-bit fixed-point representation. As can be observed from Figure [T] 
the output of the controlled system does not converge to the equilibrium point at the origin because of the 
fixed-point implementation error in the controllers. Furthermore, the practical stability region using gains K2 
and L2 is much smaller than the one using gains Ki and Li. 



Using bounds on implementation errors for the two controllers (described in Section [3]) and Proposition |2.3[ 
we can prove that the output of the system with feedback and observer gains Ki and Li (resp. K2 and L2) 
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converges to a ball centered at the origin with radius 0.5486 (resp. 0.0513), whenever the output of the system 
and the state of the observer take values in the interval [—1, 1] and the feedback gain and observer dynamic 
are implemented using 16-bit fixed-point implementation. As can be seen, given a 16-bit implementation, 
feedback and observer gains K2 and L2 may be preferred to gains Ki and Li because they have guaranteed 
bounds on practical stability region that is 10 times smaller than gains Ki and Li and provide approximately 
similar performance. If one considers the effect of disturbance and measurement noise, it can be proved that 
the output of the system with feedback and observer gains Ki and Li (resp. K2 and L2) converges to a ball 
centered at the origin with radius 5.04895(ei) -I- 0.5486 (resp. 2.53416(ei) + 0.0513), where 6(ei) is an upper 
bound on the size of the vector ei introduced in (2.6 1. 



Optimization objectives The above example suggests that the control design should optimize for the follow- 
ing objectives: the LQR and the LQG costs for performance, error caused by disturbance and measurement 
noise, and the implementation error given by a fixed-precision encoding. Accordingly, we define a cost function 
that is weighted sum of the four factors: 



(2.21) 



J{K,L)=wi 



\S{K)\\ 
\\S*\\ 



W2- 



\PiL)\\ 
\\P*\\ 



W3 



lly 



where wi, ... ,104 are weighting factors, S* and P* 



and (2.191 using standard LQR and LQG gains (Klqr and Llqg)i 
the L2 gains in (2.12) using feedback and observer gains K and L (resp 



are matrices, computed from Lyapunov equations in (2.17) 
-fly and (resp. ^{y and 73^^) are 

Klqr and Llqg) and 6(62) (resp. 

6* (62)) is the bound on the implementation error of given feedback and observer gains K and L (resp. Klqh 
and Li^Qo). Minimizing the terms and 721,^(62) inside (2.21) results in a tighter bound on the set Ay in 
Proposition 2.3 Since the four factors in (2.21 ) have different scales, we normalized them by their values using 
the standard gains K^Qfi and Llqg- The designer can choose wi, . . . , W4 based on the priorities on LQR and 
LQG performances and steady state error. Our objective is to find feedback and observer gains that minimize 
the cost function J'. 

We focus on implementation errors arising out of fixed-precision arithmetic. The bound b{e2) is computed 
using the strategy, explained in Section |3] Note that the cost function J' is not necessarily convex with respect 
to the feedback and observer gains K and L. Therefore, the proposed design strategy cannot be formulated 
as a convex optimization problem. We use a heuristic stochastic optimization approach to find feedback and 
observer gains K and L minimizing J^. 

In our exposition, we consider the plant model to be precise, and only consider quantization effects as the 
source of error. Our methodology can consider both additive and multiplicative uncertainties in the plant 
model as well [GL941. We can take those uncertainties into account by adding appropriate extra terms to the 



cost function in (2.21) using the results provided in jZSKG09l IZKGS07] . We omit the details for simplicity. 



3. Computing Quantization Error 



In this section we show how to compute a bound on the fixed-point implementation error for given feedback 
and observer gains K and L. We assume that the outputs of the controlled system and the state of the 
observer are restricted to compact subsets Y cW and X C K", respectively. 

3.1. Best fixed-point implementation. A fixed-point representation of a real number is a triple {s,n,m) 
consisting of a sign bit s G {s,u} (for .signed and unsigned), a length n £ N, and a length of the fractional part 
TO e N. The length of the integer part is n — m — 1. Intuitively, a real number is represented using n bits, of 
which m bits are used to store the fractional part. Clearly, the largest integer portion has to fit in n — m — 1 
bits. 

A variable with a fixed-point type is represented as an integer. We associate an integer variable x with the 
fixed-point representation of a real variable x. An integer variable x that represents a fixed-point variable with 
type (u, n,TO) can be interpreted as the rational number 2~™x. We deal with a signed number by separately 
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tracking the sign and the magnitude, performing the operations on the magnitudes using unsigned arithmetic, 
and finaUy putting the appropriate sign bits back. 

An operation using real arithmetic may have different fixed-point implementation depending on how many 
bits are allocated to hold the integer part and the fraction part of the variables. Allocating fewer number of 
bits than required to hold the integer part may lead to overflow. On the other hand, if more than the required 
number of bits are allocated to the integer part, the quantization error increases due to assigning few bits to 
the fractional part. When we compare fixed-point implementation of different controllers, we first synthesize 
the best possible implementation of a controller. 

Let us fix the number of bits to be n for the implementation of a controller. With n bits, let h be the 
upper bound on the quantization error in a fixed-point implementation / of a controller for a given range 
for the inputs. The fixed-point implementation / is the best implementation if there does not exist another 
implementation /' using n bits, for which the upper bound on the quantization error is b' and b' < b. 

If the ranges of the variables in the real arithmetic computation can be computed exactly, it is possible 
to synthesize the best fixed-point implementation. In the best fixed-point implementation, the number of 
bits allocated to the integer part is just enough to hold the integer part of any value in that range. For 
example, if the range of a variable is [-35.55, 48.72], the datatype for the variable in the best 16-bit fixed-point 
representation is (1, 16, 9). 

The range computation problem of variable y in an operation y = /(xi, . . . , a;„) involves solving a maximization 
and a minimization problem where / is the objective function and the ranges on xi, . . . ,x„ form the set of 
constraints. If the function / is convex, the range of y can be computed exactly, and also it is straight-forward 
to find the best fixed-point implementation for the operation. 



3.2. Error bound computation. We apply mixed-integer linear programming based optimization technique 
to find out the error bound between a computation in real arithmetic and its best fixed-point implementation. 
Suppose we have an arithmetic operation s : a = b op c, where op € — , *}. If op = *, then either 6 or c is 
a constant. If op = -I- or op = — , then b and c can both be variables. We associate an integer variable x with 
the fixed-point representation of a real variable x. Let the range of the values for a and b and c are [Za,Wa], 
[/b,Wf,], and [lc,Uc], respectively. Let the fixed-point representation of a, b and c are {s,na,ma) , (s,n(,,TOh), 
and (s. Tic, nic), respectively. Let 6(eb) and b{ec) be bounds on the quantization errors of b and c, respectively. 
The optimization problem to find the bound on the error is given by: 



maximize 
Subject to 



(3.1) 



|a- 2-'"»a 

la < a<Ua 

Ih < b < Uf, 
h _ 2-™" * b 
c - 2-'"- * c 

a = b op c 

$(fp(5)) 



< b{eb) 

< He,) 



where fp(s) is the fixed-point representation of the statement s and $(s) denotes a logical formula that relates 



the inputs and outputs of the fixed-point representation s. Technically, $ is the strongest postcondition Wi n93j 
of s with respect to true. We compute $ using an arithmetic encoding of a fixed-point computation [AMSTlOj . 
Here we illustrate the computation of the strongest postcondition $ using an example. 

Example. Suppose we have the following arithmetic operation 

s : y = -7.2479 * x . 

Assume the compact set for x is [-1, 1]. The fixed-point expression corresponding to s in the best fixed-point 
implementation is 

fp(s) : -y = (-115 * > 6 . 
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The strongest postcondition <f>(fp(s)) of fp(s) is given by: 

$(fp(s)) := tmp = -115* X A 

trap > — >■ tmpl = trap A 

trap < — > trapl = —trap A 

tmpl = 2^ * divisor + remainder A 

remainder > A remainder < 2^ A 

tmp > ^ i) — divisor A 

tmp < ^ y — —divisor , 

where tmp, tmpl, divisor, and remainder are integer variables. 

Depending on the arithmetic operation, we need to solve at most 4 instances of mixed integer linear pro- 



gramming problem to solve the optimization problem in (3.1 ), and the maximum among all of them gives the 
bound on the error in the fixed-point implementation. 

We use the above technique to compute the bound on the error in one operation in the fixed-point implemen- 
tation of a gain. The implementation of a gain involves a series of arithmetic operations. We compute the 
error bound for the output of one arithmetic operation at a time. Let s : a = 6 op c is an arithmetic operation 
in the implementation of a gain. In the arithmetic operation, b and c may either be a constant, a state 
variable or a temporary variable which captures the result of some previous operation. If b (or c) represents a 
constant, and the fixed-point representation contains m bits for the fraction part, then the error in the fixed 
point representation is bounded hy If b (or c) represents a state variable, then the fixed-point datatype 
can be determined from the given compact set for the state, and the fixed-point datatype can be determined 
accordingly. Then the error in the fixed-point representation is bounded by ^ , where m is the number of bits 
to represent the fraction part in the fixed-point datatype of the variable. If b (or c) is a temporary variable 
used to hold the result of an earlier computation, then the range and error bound for the variable is already 
known. 

4. Optimal Controller Synthesis 



We now describe our controller synthesis algorithm that minimizes the cost function (2.21) combining LQR 
and LQG performance, disturbance, measurement noise and implementation errors. Since the cost function is 
non-convex, we use a stochastic local search technique. 

4.1. Particle swarm optimization (PSO). PSO is a stochastic local search approach. It maintains a set of 
potential solutions (called "particles") in a compact d dimensional search space D = 0^=1 bmin' J/max] C M.'^, 
minimizing a given cost function. The particles move in this space according to their velocity. Each particle, 
indexed by i G N, has a position i/i E M'^, changing between j/min and j/max, and a velocity vector Vi € M**, 
changing between some vectors Umin and Wmax- The terms Vmin and Wmax are often set to the maximum 
dynamic range of the variables on each dimension jZKGSPOO] : — = t^4ax = l2/max ~ ymini- Every particle 
remembers its own best position (i.e., the lowest value of the cost function achieved so far by this particle) 
in a vector Pi G M'^. The best position with respect to the cost function among all of the particles so far is 
stored in a vector Pg gM.'^. 

PSO updates the positions and velocities of all particles iteratively. The new velocity and position for particle 
i are determined as: 

(4.1) vl+' ^w'v\ + cin {Pi - yl) + C2r2 (P^ - 2/^ ' 

(4.2) y^i=y^+.^\ 

where the superscript I denotes the iteration number, the subscript i = 1,. . . ,N denotes the index of the 



particle, and N is the number of particles. The constant w in (4.11 is updated using the inertia weight 
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approach [EKG12) as the following: 

(4.3) w = max <^ w,„in, w^ax , (« - 1) 

where Wmax and Wmin are adjusted to 1 and ^^^^ — 1 and /max is the maximum number of iterations. The 



constants ci and C2 in (4.1 1 are the acceleration constants, influencing the convergence speed of particles 



toward its own and global best positions and set to 0.5 and 1, respectively |EKG12] . The constants ri and r2 



in (4.1 1 are uniformly distributed random numbers on the interval [0, 1]. 



4.2. Overall algorithm. The PSO algorithm is used to search for feedback and observer gains K e 



and L g W^^p for the control system (2.5), minimizing (2.21 1. Note that a particle in PSO represents a 
feedback and an observer gain K and L, respectively, moving inauTOXn + nxp dimensional search space. 
To discard those gains that make the controlled system unstable, we penalize unstable gains by including a 
penalty term P in the cost function such that P = Q \i — B^-K and — LC are Hurwitz and P — +oo 
otherwise. The cost function for PSO is then F{K, V) = J{K, L) + P{K, L). 

The design steps can be summarized as the following: 

(1) Initialize positions of N feedback gains Ki and observer gains Li by Klqji and Llqg, respectively, 
and uniformly randomly initialize their velocities , i = 1, . . . , iV, where N is the number of particles. 

(2) Given any initial feedback gain Ki and observer gain Li, compute the cost function F{Ki, Li). To 
compute P, check if Ar — BrK and Ar — LC are Hurwitz. There are some steps to compute J . First 



compute S{Ki) and P{Li) by solving the Lyapunov equations (2.171 and (2.19), respectively, and find 
their induced 2 norm. Second, compute the £2 gains 7iy and 72^. Third, compute 6(62) by solving 
the optimization problems from Section [3j 
(3) Compare F{Ki, Li) to its own best position Pi so far and the global best position Pg so far. If F{Ki,Li) 
is less than the previous personal best (resp. the global best), update the best position (resp. the global 
best) to Ki and Li. 



(4) Modify the velocity and position of each pair Ki and Li according to (4.1) and (4.2 1. 

(5) If the number of iterations, denoted by I, reaches the maximum, denoted by /max, or the value of F 
does not change for the global best position Pg for 50 consecutive iterations up to error 10~^ then go 
to Step (6), otherwise go to Step (2); 

(6) The latest Pg is the optimal controller. 

5. Extension: PID Controllers 

FID controllers are a common class of controllers in many industries, such as automotive, power systems, ser- 
vomotors, and so on. We now extend the analysis of Section[2]to PID controllers. A PID controller generalizes 
a proportional feedback controller, and includes three terms: proportional, integrator, and differentiator. For 
an input u, the output 1] of the PID controller is computed as follows: 



(5.1) ri{t)^Kpv{t)+Ki I v{s)ds + KD—j^,yteR^, 



* dv{t) 



where Kp, Kj, and Ko are called proportional, integrator, and differentiator gains, respectively. To describe 
the mismatch between the PID specifications and its software implementation, we consider the discrete-time 



version of (5.1). A common way of discretizing an integrator is based on the trapezoidal approximation. An 



integrator term: 

r]{t) = I v{s)ds, yt e M(j", 

can be discretized as follows: 



T 



(5.2) y[r + 1] = y[r] + - {u[r + 1] + u[r]) , Vr G No, 
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where r is the samphng time, y[r] — r](rT)-\-ei and u[r] — v^rr), for any r e Nq. A common way of discretizing 
a differentiator, is based on the backward Euler method. A differentiator term: 



dv{t) 
dt '■ 



can be discretized as foUows: 
(5.3) 



u\r + 11 — u\r] , , 
y[r + 1] = ^ J Li^ Vr e No, 



where y[r] — ri{rT) + 62 and u[r] — u(rr), for any r G Nq. By using fast samphng time assumption, we 
can ignore the errors ei and 62 in the discretized versions of integrator and differentiator in comparison with 
quantization errors. To follow the same analysis as in Section [2j we need a state space realization of PID 
controller. By resorting to results in control classic [KaiSOj and using the discretization rules in (5.2 1 and 
(5.3), the state space realization of discretized PID controller with input u[r] and output y[r] is obtained as 



follows: 

(5.4) 

where 



x[r + 1] = Ax[r] + Bu[r], 
y[r] — Cx[r] + Du[r], 



" 1 ■ 




■ ■ 




, B = 




1 


1 



Kd ^ Kd 

KjT 

r T 



D = K 



Kit Kd 



Without loss of generality, consider a single-input {m — 1) single-output (p = 1) discrete-time linear control 
system of the form: 

x[r + 1] = Ax[r] + Bu[r], 
y[r] = Cx[r]. 

Since the input of the PID controller is equal to the negative of the output of the plant {u = —y) because of 
negative feedback and the output of the PID controller is equal to the input of the plant {u — y), one obtains: 



(5.5) 



x[r + l] = \^A- BDCj x[r] + BCx[r], 
x[r + 1] = -BCx[r] + Ax[r]. 



Similar to what explained in Section [2] by fixed-point implementation of the PID controller, one gets the 
following overall dynamic: 



(5.6) 



x[r + 1] = {^A- BDCj x[r] + BCx[r] 
x[r + 1] = -BCxlr] + Ax[r] + e^i, 



Be 



q2, 



where e^i and 6^2 are quantization errors in computing the PID controller. Now, we can use the same 
strategy, as explained in Subscction |4.2[ to design parameters Xp, Ki, and Kd of PID controllers minimizing 
a performance-based cost function as well as the effect of quantization error. For example, one can consider: 



PM 



W2 

GM 



- W37(6(eqi) +6(6,2)), 



(5.7) J{Kp,Ki,Kd) 

where PM and GM are phase and gain margins, Wi,W2, are weighting factors, 7 is the £2 gain of the linear 
control system (5.6) and b{eqi) and 6(6,2) are the bounds on the implementation errors 6,1 and 6,2- Note that 



control over PM and GM guarantees robust stability of the closed-loop systems |Hes09| . The phase and gain 
margins measure the system's tolerance to the time delay and the steady state gain, respectively. 
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Control systems 


# bits 


Synthesized gains 


Time cost 


K 


L 


Bicycle 


16 


[3.0253 12.6089] 


[0.0132 0.1021]" 


Ih36m41s 


DC motor position 


16 


[0.1129 0.0211 0.0093] 


[0.0390 0.3700 - 0.0175]-' 


Ih39m06s 


Pitch angle control 


32 


[ -0.1202 42.5655 1.0001] 


[0.0001 0.0017]-" 


8h31m53s 


Inverted pendulum 


32 


[-1.5362 -2.0254 16.5192 2.7358] 


r 0.0017 0.0021 0.0012 l ' 
0.0001 0.0018 0.0122 0.0770 


9h54ml7s 


Batch reactor process 


16 


r 0.0583 0.9093 0.3258 0.8721 ' 
-2.4638 -0.0504 -1.7099 1.1653 


r 0.0774 -0.0022 0.0267 0.0356 1 
-0.0103 0.0227 0.0398 0.0001 


3h08m29s 



Table 1. Synthesized gains and required time for synthesizing them. 



Control 
systems 


lub of LQR cost 


LQG cost 


Steady state error 


LQR 


Synthesized 
K 


LQG 


Synthesized 
L 


LQR-LQG 


Synthesized 
gains 


Bicycle 


3956.3||a;||^ 


4331.7||x||^ 


0.0229 


0.0246 


5.04896(ei)-|-0.5486 


2.53416(ei)-|-0.0513 


DC motor position 


1001.6 a; 


1376. 7 ^^" 


36.6315 


36.6731 


30.5666(ei)-|-0.16 


15.42i;)(ei)-f0.011 


Pitch angle control 


2.9732 X lO^llxll^ 


2.9887 X lO^llxll^ 


0.0013 


0.0018 


2.67816(ei)-|-0.4746 


1.44536(ei)-|-0.0807 


Inverted pendulum 


4.2988 X 10* x ^ 


5.3471 X 10* X 


0.3600 


0.3897 


83.4217fe(ei)+0.0432 


30.380i;)(ei)+0.0086 


Batch reactor process 


223.1773||x||^ 


223.1825||x||^ 


0.0731 


0.0949 


2.93096(ei)-|-0.4194 


2.12166(ei)-|-0.1642 


Table 2. Least upper bound (lub) on the 


LQR cost (2.131, for a given initial condition x, 



the LQG cost (2.14), and the Euchdean norm of the steady state error for the LQR-LQG and 



the synthesized gains. 



6. Experimental Results 



We implemented the algorithm presented in Section 4.2 in Matlab. We use a PSO function in Matlab, 
introduced in [EKG12j . We implemented a static analyzer in Ocaml that synthesizes the best fixed-point 
program and computes the bound on the fixed-point implementation error for given feedback and observer 
gains K and L, respectively. The tool gets the number of the bits in the fixed-point datatype, compact 
subsets y C and X C K", and feedback and observer gains K and L, respectively, as inputs. The 
optimization problems in computing the error bound are solved using the mixed-integer linear programming 
tool lp_solve 



We applied the proposed controller synthesis approach to a number of linear control systems. All the experi- 
ments were done on a laptop with CPU Intel Gore 2 Duo at 2.4 GHz. In all of the experiments, the number 
of the particles in PSO is iV = 24, the m aximum number of iterations is se t to Zmax — 100, and we choose the 
matrices Q = /„, and R — 1^ in (2.13) and Q — Iq, and R = Ip in (2.16). The value of Zmax was chosen in 



such a way that appropriate gains are obtained in terms of the cost function (2.21) (or (5.7)) for all control 



systems. JVIoreover, we assume that the search space is D — J|"^-|™ [— 150, 150] C iu"><™+"><p that is large 
enough and contains the standard LQR and LQG gains for all the examples. Furthermore, without loss of 
generality, we work on the compact subsets Y — YVi=i[^^' 1] C K'' and X = n"=i[~l' 1] ^ K". AH constants 
and variables are expressed in SI units. 



Our unstable examples include a bicycle [AM08| . a DG motor position control jcilltij . a pitch angle con- 
trol jCmuj ■ an inverted pendulum [Cmuj . a batch reactor process |GL94j and another inverted pendulum for 
PID synthesis [Cmtlj . See Table [l] and [2] for experimental results. Note that for those examples for which 32-bit 
implementation is chosen, the 16-bit one provides a stability region which is even larger than the range of the 
variables inside the controller. As can be seen from Table [2l in comparison with the conventional LQR-LQG 
approach, the proposed synthesis approach in this paper worsens the LQR and LQG performances by at 
most 1.37 times (for DG motor position) and 1.38 times (for Pitch angle control), respectively. However, the 
proposed synthesis approach improves the size of the region of practical stability due to quantization error by 
at least 2.55 times. For certain examples, the improvement goes beyond the factor of 10. For bicycle and DG 
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Best: 3.1406, Mean: 3.1605 



o 











Best 
-Mean 





















































20 40 60 80 100 

Iteration 



Figure 2. Cost of the best particle and average cost of all population vs iteration. 



motor position, the region of practical stability due to quantization error improves by a factor of 10.69 and 
14.55, respectively. 



The detailed description of the systems are available as follows. 



W2 



Bicycle The model of a bicycle is shown in (2.20). The weighting factors in (2.21) are chosen as Wi 

5. The results of the LQR, LQG and the proposed method are shown in Tables [T] and |2] 



W3 = 1 and W4 

Figure [2] shows how the value of the cost function improves with the number of iteration. The figure shows 
how the value of the cost function monotonically decreases with the number of iterations. The fixed-point C 
code for the synthesized controller is shown in Figure [3j 

DC motor position control The dynamic of a DC motor position control, borrowed from fcmui, is given 
by: 
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= [1 0] 
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K 






Jk 
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where is the angle of the motor's shaft, ^2 is the angular velocity of the motor's shaft, ^3 is the armature 
current, h = 3.508 x 10^^ is the damping ratio of the mechanical system, J = 3.228 x 10~^ is the moment 
of inertia of the rotor, K — 0.027 is the electromotive force constant, i? = 4 is the el ectric resistance, 
L = 2.75 X 10~^ is the electric inductance and v is the source voltage. The weighting factors in ( 2.21 ) are chosen 
as wi = W2 ~ = 1 and W4 — 5. The LQR and LQG gains are given by Klqr = [0.4055 0.3782 0.0022] 
and Llqg = [0.0288 0.3858 — 0.0026]-^ and the gains, computed by the proposed approach in this paper, are 
given in Table [T] The detailed results are shown in Tables [T] and [2] 



Pitch control The dynamic of the longitudinal motion of an aircraft, borrowed from [emu], is given by: 

(V+LO) 







-0.313 


56.7 






-0.0139 


-0.426 


. 4 . 









56.7 








" a " 








1] 












6 







" 6 " 




0.232 




6 


+ 


0.0203 




6 








where ^ is the angle of attack, ^2 is the pitch rate, ^3 is the pitch angle, and v is elevator deflection angle. The 



weighting factors in (2.21 ) are chosen as wi = W2 ^ W3 = 1 and — 5. The LQR and LQG gains are given 
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float output (float yin) 
{ 

static int xl = x\q; 1 1 f ixdt (1 , 16 , 14) 

static int x2 = x2o; II f ixdt (1 , 16 , 14) 

int xl.neiu; // f ixdt (1 , 16 , 14) 

int xljnew; II f ixdt (1 , 16 , 14) 

int ti; // fixdt(l,16,ll) 



// Intermediate variables 

int Gain\\ II f ixdt (1 , 16 , 15) 

int Gaini; II f ixdt (1 , 16 , 15) 

int Gaini; II f ixdt (1 , 16 , 15) 

int Add\\ II f ixdt (1 , 16 , 14) 

int Gom4; // f ixdt (1 , 16 , 15) 

int Gairia; II f ixdt (1 , 16 , 15) 

int Gom6; // f ixdt (1 , 16 , 15) 

int Addl; II f ixdt (1 , 16 , 15) 

int Gain7; II f ixdt (1 , 16 , 13) 

int GainS; II f ixdt (1 , 16 , 11) 



y = convert_to_f ixedpoint (ym) ; 

Gaini = (31499 * xl) » 14; 

Gain2 = (-3145 * a;2) >> 14; 

Addl = (Gaini + Gain2) » 1; 

Gaini = (432 * y) >> 14; 

xl_new = ({Addl << 1) + Gaini) >> 1; 

Gaini = (-1907 *xl) >> 14; 

Gain5 = (23835 * x2) » 14; 

Add2 = Gaini + Gain5; 

GainQ = (3345 * y) » 14; 

x2_new = (Addl + GainQ) » 1; 

Gain! = (24783* xl.new) >> 14; 

Gains = (25823 * x2_new) >> 14; 

u = (Gain! + (GainS << 2)) » 2; 

returnCf loat Cu) ) ; 



Figure 3. synthesized fixed-point controller C code for Bicycle. 



by Klqr = [-0.1141 49.1428 0.9995] and Llqg = 10"^ x [0.6407 0.0039 O.eeSS]"^ and the gains, computed 
by the proposed approach in this paper, are given in Table [T] The detailed results are shown in Tables [l] and 

HI 



Inverted pendulum Consider a simple physical model of an inverted pendulum on a cart, borrowed from 
[cniuj . The dynamics of the system is given by: 
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where ^i, and ^2 are the position and velocity of the cart, respectively, ^3, and ^4 are the angular position 
and velocity of the mass to be balanced, v is the applied force to the cart, g = 9.8 is the acceleration due to 
gravity, / = 0.3 is the length of the rod, m — 0.2 is the mass of the system to be balanced, M = 0.5 is the 
mass of the cart, 6 = 0.1 is the coefficient of friction of the cart, and / = 0.006 is the inertia of the pendulum. 
The weighting factors in (2.21) are chosen as wi = ii;2 = W3 = 1 and W4 = 5. The LQR and LQG gains are 
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given by Klqr = [-0.9929 - 2.0276 20.2819 3.9126] and 



L 



LQG 



0.0016 0.0011 0.0007 0.0034 
0.0007 0.0051 0.0111 0.0618 



and the gains, computed by the proposed approach in this paper, are given in Table [T] The detailed results 
are shown in Tables [l] and [21 



Batch reactor process Consider an unstable batch reactor process, borrowed from [ GL94j . The dynamic of 
the system is given by: 

^1 

^4 



1.38 


-0.2077 


6.715 


-5.676 
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-0.5814 


-4.29 





0.675 






1.067 


4.273 


-6.654 


5.893 






0.048 


4.273 


1.343 


-2.104 




$4 











" 1 " 


5.679 





V + 
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1.136 


-3.146 


1 


1.136 







1 



■n = 



1 
1 



^4 



The weighting factors in (2.21) are chosen as wi 
are given by: 



W3 = 1^ W2 — 2, and W4 = 5. The LQR and LQG gains 



K 



LQR 



^LQG 



0.0376 0.9157 0.3262 
-2.4884 -0.0734 -1.7461 



0.0447 -0.0003 0.0170 0.0127 
0.0020 0.0058 0.0059 



0.8226 
1.1438 

T 



and the gains, computed by the proposed approach in this paper, are given in Table [T] The detailed results 
are shown in Tables [T] and [21 

PID controller In this example, we provide a PID controller for an inverted pendulum whose dynamic is 
given by a transfer function. Consider the transfer function of an inverted pendulum, borrowed from jcmuj . 
given by: 



(6.1) 



Ms) 



U{s) 



q 



{A'l-\-m)mgl 



bmgl 

q " q 

where q = {M + m){I + mP) — (ml)'^, output is the angular position of the mass to be balanced, input v 
is the applied force to the cart, g = 9.8 is the acceleration due to gravity, I = 0.3 is the length of the rod, 
TO = 0.2 is the mass of the system to be balanced, M = 0.5 is the mass of the cart, b = 0.1 is the coefficient 
of friction of the cart, and / = 0.006 is the inertia of the pendulum. Using standard results in control theory 
[Kai80j . one obtains the following state space realization for the inverted pendulum: 
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Our objective is to design PID gains Kp, Kj, and Kj) minimizing the cost function (5.71 with weighting factors 
Wi — W2 = = 1 and the closed loop system has a settling time (tg) of less than 5 seconds and pendulum 
should not move more than 0.05 radians away from the vertical axis. The latter two constrains are treated the 



same as the stability constraint in Subsection 4.2 by penalizing the cost function (5.7). The synthesized gains 



are Kp = 109.032, Ki = 1.2268, and Kd = 13.9945. The closed loop system has PM +00, GM = 26237, 
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7(6(eqi) + b{eq2)) = 4.1705 x 10 ^, settling time tg = 0.4790, and pendulum does not move more than 0.0098 
radians away from the vertical axis. 



7. Conclusion 

We have presented a generic methodology to search for optimal controller implementations that minimize 
implementation errors in addition to traditional controller performance criteria. While we have instantiated 
the methodology using the LQR and LQG costs and quantization errors, our algorithm is more generally 
applicable to other performance criteria and other sources of modeling or implementation error. The controller 
synthesis algorithm can be seamlessly added to any design and automatic code generation tool flow to enhance 
its capability to generate correct-by-construction high performance controller software. By automatically 
synthesizing the minimal error controller, we sidestep the need for post-design verification. 
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